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ABSTRACT 

We study the chaotic orbital evolution of planetary systems, focusing on secular (i.e., orbit- averaged) 
interactions, because these often dominate on long timescales. We first focus on the evolution of a 
test particle that is forced by multiple massive planets. To linear order in eccentricity and inclination, 
its orbit precesses with constant frequencies. But nonlinearities modify the frequencies, and can shift 
them into and out of secular resonance with the planets' eigenfrequencies, or with linear combinations 
of those frequencies. The overlap of these nonlinear secular resonances drive secular chaos in planetary 
systems. We quantify the resulting dynamics for the first time by calculating the locations and widths 
of nonlinear secular resonances. When results from both analytical calculations and numerical inte- 
grations are displayed together in a newly developed map, the "map of the mean momenta" (MMM), 
the agreement is excellent. This map is particularly revealing for non-coplanar planetary systems and 
demonstrates graphically that chaos emerges from overlapping secular resonances. We then apply this 
newfound understanding to Mercury. Previous numerical simulations have established that Mercury's 
orbit is chaotic, and that Mercury might even collide with Venus or the Sun. Guided by intuition 
from the test particle case, we show that Mercury's chaos is primarily caused by the overlap between 
resonances that are nonlinear combinations of four modes, the Jupiter-dominated eccentricity mode, 
the Venus-dominated inclination mode and Mercury's free eccentricity and inclination. Numerical 
integration of the Solar system indeed confirms that a slew of these resonant angles alternately librate 
and circulate. We are able to calculate the threshold for Mercury to become chaotic: Jupiter and 
Venus must have eccentricity and inclination of a few percent. Mercury appears to be perched on the 
threshold for chaos. 



1. INTRODUCTION 

The question of the stability of planetary orbits in the 
Solar system has a long history, and has attracted the at- 
tention of some of the greatest scientists, including New- 
ton, Laplace, Lagrange, Gauss, Poincare, Kolmogorov, 
and Arnol'd. Newton thought that interplanetary per- 
turbations are eventually destabilizing, and that divine 
intervention is required to restore the planets' orbits to 
their rightful places (Laskar 1996). Yet it is only over the 
last twenty years that the stability of the Solar system 
has been def initively settled, with the aid of computer 
simulations (Sussman fc Wisdom 1988 Laskar 1989| 



I Quinn et aL]rl991t [Wisdom & Holman|ri99l[ |Lecar et al.] 
|2001t ILaskar & Gastineau||2009p. We now know that 
Newton was not far oft": the Solar system is marginally 
stable: it is unstable, but on a timescale comparable to 
its age. In the inner Solar system, the planets' eccentric- 
ities chaotically diffuse on a billion- year timescale, with 
the two lightest planets. Mercury and Mars, experienc- 
ing particularly large variations. In fact. Mercury has 
roughly a 1% chance of colliding w ith Venus or the Sun 



withi n the next five billion years (Laskar & Gastineau 
2009 ) . By comparison, the giant planets in the outer So 



lar system are well-spaced, and their orbital elements un 
dergo largely quasiperiodic variations, exhi biting chaotic 
diffusion only on extremely long timescales ( Laskar|1996 
Murray fc Holman 1999). 



'I'he chaotic dynamics in the inner Solar system is pri- 
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marily due to ^eci^kr interactions ( Laskar|2008 ) . In gen- 
eral, interplanetary interactions can be decomposed into 
secular ones and MMR's (mean motion resonances — not 
to be confused with secular resonances). Secular interac- 
tions result from orbit-averaging the equations of motion. 
Since averaging a Keplerian orbit produces an elliptical 
ring, secular evolution can be thought of as interactions 
between elliptical rings. Secular timescales are long — 
they are longer than the orbital time by at least the ra- 
tio of the star's mass to that of a planet. By contrast, 
interactions driven by MMR's depend on orbital phase, 
and typically occur on the orbital timescale, or longer if 
some of the planets' orbital periods are close to integer 
ratios. Intuitively, one would expect that the dynamics 
on long timescales can be treated by averaging over the 
fast orbital phase — i.e., that they are secular in nature. 
This is true in the inner Solar system. It is also true more 
generally for well-spaced planets that do not happen to 
lie near mean motion resonances 

Linear secular theory has been understood for hun- 
dreds of years, dating ba ck to the famous solu tion of 
Laplace and Lagrange (see 'Murray fc Dermott"2000). To 
linear order in the planets' eccentricities and inclinations, 
secular theory reduces to a simple eigenvalue problem, 
with two eigenmodes per planet — one for the eccentric- 
ity degree of freedom, and one for the inclination. Each 
eigenmode has a constant amplitude and a longitude that 
precesses uniformly in time. But linear secular theory is 
clearly incapable of describing the chaotic orbits of the 

^ In the outer Solar system the dynamics is not mainly secular 
because the giant planets lie near a number of MMR's, such as the 
5:2 between Jupiter and Saturn (the "Great Inequality"), and the 
2:1 between Uranus and Neptune. 
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Solar system. 

Despite the importance of secular chaos, there has 
been surprisingly little theoretical understandin g of it 
(see, e.g., the review of Solar system chaos by Lecar 
By contrast, chaos due to MMR's is well- 



let al. 2001) 

understood, a nd accounts fo r the Kirkwood gaps in the 
asteroid belt ( Wisdom] 1983 ), and for the very weak chaos 
of the ou ter Solar system planets, which is due to 3-body 
MMR's ( [Murray fc Holman|1999| ). In ah cases that have 
been studied in the Solar s ystem, chaos is caused by over- 
lapping resonances (Chiri kov|1979[ [Lecar et"ar|2 001). In 
linear secular theory, there can be secular resonances. 
And it is generally supposed that the chaos in the in- 
ner Solar system is caused by overlapping secular reso- 
nances. Yet thus far there has been little quantitative 
calculation. To our knowledge, the only previous the- 
oretical work towards calculating secular chaos was by 



Sidlichovsky (1990) , who considered the coplanar case, 
as we describe below (Q. 
Numerical attempts to identify the mech anism o f chaos 



1992) and Sussman & Wisdom 



in th e inn er Solar system were made by [Laskarj (1990, 

( [T992I ). 'I'hese autliors 
tound that the angle associated with the (secular) fre- 



quency (^mercury ^jupiterj 



mercury 



g) alter- 



nately librated and circulated in their simulations, where 
g is the apsidal precession rate, and s is the nodal pre- 
cession rate (or, to be more precise, g and s here re- 
fer to the frequencies of the normal mode that is domi- 
nated by the corresponding planet). Laskarj ( |1992|) also 
found that two angles associated with Earth and Mars, 

corresponding to 2(^mars - dearth) - (Smars " ^earth) and 

(^mars - dearth) - (^mars - ^earth) , alternately librated, and 
conjectured that the overlap of tho se secular resonances 
was responsible for chaos. But, as [Sussman fc Wisdom 
( [1992 ) note, Laskar's conjecture is not fully convincing, 
because there are too many unrelated angles that alter- 
nately circulate and librate, and it is not clear which are 
dynamically important. Furthermore, only one librating 
angle has been identified for Mercury. Yet chaos requires 
the over lap of a t least two resonances, so why is Mercury 
chaotic ( Lecar et al.|20 01)? Without a theory for secular 
chaos, the dynamics remain obscure. For example, why 
does instability in the Solar system occur at such low 
values of eccentricity and inclination (~ few percent)? 
What sets the timescale of the chaos? Can secula r chaos 
shape the architecture of the inner Solar system (Laskarj 
1^96)? And can it shape the architecture of extrasolar 



planetary systems (Wu fc Lithwick 2010)? Without a 



theory for secular chaos, we will be forever at the mercy 
of computer simulations. 

In this paper, we construct the theory for secular chaos 
of a test particle, and then apply the theory to Mercury. 
In pi we present the test particle's equations of motion. 
In |3[ we describe the coplanar solution, and in Q we 
generalize to the case when bodies have non-zero inclina- 
tions. In ^ we apply the theory to N-body simulations 
of the real Mercury. We conclude in ^ 

2. SECULAR EQUATIONS OF MOTION 

We focus on the secular evolution of a massless test 
particle that is orbiting a star in the presence of multiple 
massive planets, assuming the planets' orbits are known. 

The particle has six orbital element s, |a, e, i, A, 
using standard notation (Murray fc Dermott 2000[ ). In 



secular theory, one averages over A. As a consequence, 
a is a constant of motion, leaving only four orbital 
elements to be considered. The equations of motion 
for the particle's eccentricity and longitude of periapse 
(e and m) are given by Hamilton's equations for the 
Poincare cano nical variables F = y /GMQa (l — \/r^-^) 
and 7 = —w (Murray fc Dermott 2000). Since a is con- 



stant, it is simpler to choose the canonical momentum to 
be oc T / ^GM(^a^ so we introduce the momentum 



= 2(1- Vl-e') 
:e2 + (9(e^) 



(1) 

(2) 

and its conjugate co-ordinate to be w. Although this is a 
non-canonical transformation from Poincare's variables, 
if we simultaneously re-scale the energy by defining as 
the Hamiltonian 



H 



(3) 



where E is the particle's energy per unit mas^then the 
equations of motion are Hamilton's equations. 



dvu/dt = dH/dpe 
dpe/dt = —dH/dw 



(4) 
(5) 



Therefore we may consider (pe^ru) to be canonically con- 
jugate. Similarly, for the inclination and longitude of 
node (i, 1^), we take the canonical variables to be related 
to the corresponding Poincare variables in the same way 
by defining 



Pi=Wl-e'^ sin^(z/2) 



(6) 



and taking its conjugate co-ordinate to be 1^. The equa- 
tions of motion for {pe^vu^pi^ Q) are Hamilton's equations 
generated by the scaled H. These equations are exact as 
long as a is constant, which is the case for secular inter- 
actions. 

An alternative formulation of the equations of mo- 
tion will also prove useful. For an arbitrary Hamilto- 
nian H{p, g), one may define the complex canonical vari- 
ably Z = y^e*^. As may be directly verified, the equa- 
tion of motion for Z is then dZ/dt = idH/dZ"^ where 
H{Z^Z^) = H{p^q). This complex equation of motion 
simultaneously encodes both of Hamilton's equations. 
Since our real canonical variables are (pe^ru^Pi^ft), we 
introduce the complex ones 



/Pee 



2 1 



/Pie ' 



in 



2(1 



2\l/4 



1/2 



sin(i/2)e 



^ ee 



(8) 
(9) 



which, to leading order in e and i, are the usual complex 
eccentricity and inclinationj^ The equations of motion 

^ The test particle's energy per unit mass E is given by equation 
([All) in Appendix A for the case of a single external planet. 

Our defin ition of the complex canonical variable differs from 



lOgilviel (12007) by a minus sign, and hence our Hamilton's equation 
also ditters by a minus sign. 

^ The symbol e denotes both the eccentricity and the exponen- 
tial (Euler's constant), and the symbol i denotes both the inclina- 
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for {z, Q are 



dz/dt = idH/dz* 
dC/dt = idH/dC 



(10) 
(11) 



Throughout this paper, we freely switch between the set 
of real canonical variables and the set of complex ones, 



RCV : {pe.Tu^pi.Q) 
CCV: (^,0 



(12) 
(13) 



Although Hamilton's equations of motion for the RCV 
and CCV are exact, we take an approximate form for 
the Hamiltonian by expanding H to fourth order in e 
and z, and to leading order in the ratios of semimajor 
axes, keeping only secular terms. The relevant terms are 
listed in Table [T] in Appendix A. Throughout the bulk of 
this paper, we focus on the ci, C2, C4, cn, C 12, C14, and C15 
terms in that table. As we show in §4.5[ the remaining 
terms in the table have a small effect in the parameter 
regime we focus on. 

3. COPLANAR JUPITER AND SATURN 

In this section, we consider the evolution of a test 
particle in the presence of two massive exterior planets 
("Jupiter" and "Saturn"), with all bodies having zero in- 
clinations. We assume that the massive planets' orbits 
are given by their linear Laplace-Lagrange solution, and 
evolve the test particle's equations of motion to leading 
nonlinear order (terms listed i n Table [ l ] in A ppendix [A|). 
This case was worked out by Sidlichovsky (i990). We 
describe it here in some detail because it sets the stage 
for the considerably more complicated case with non-zero 
inclinations (Q. 

3.1. Jupiter Only 

We shall solve the coplanar case with a sequence of 
increasingly complicated sub-cases. Consider first the 
linear secular evolution of a test particle perturbed by 
a circular Jupiter. From Table [l] in the Appendix, the 
particle's Hamiltonian is 



H{z)=^\z\ 



where the constant 



7 ^ 



3 mj 3 / 
a 

4 Mo V 



1/2 



(14) 



(15) 



is the particl e's l inear free precession rate induced by 
Jupiter (eq. A4 ); mj is Jupiter's mass, and a is the 
ratio of the particle's semimajor axis to Jupiter's. The 
equation of motion is dz/dt = idH/dz^ = ijz^ with solu- 
tion z = const X e*^^, an orbit with constant eccentricity 
that precesses at frequency 7. 

For the second sub-case, we consider the test particle's 
linear evolution when Jupiter is assigned a constant ec- 
centricity ej and precession rate ^j, in which case the 
particle's Hamiltonian is 



H{z)=^{\2 



(eje 



igjt . 



■c.c.)) 



(16) 



tion and the imaginary unit. There should be no confusion because 
for the remainder of this pape r we use as our dynamical variables 
either the RCV or CCV (eqs. [Tsl-fTa]) in lieu of e and i. 



where ej = \aej] we drop nonlinear terms (i.e. the 
fourth-order terms in Table [T]). Of course, if Jupiter were 
the only massive planet in the system, it would not pre- 
cess {gj = 0). But we consider a finite gj in anticipation 
of what happens when Saturn is included. 
Hamilton's equation is 



1 dz 
ij dt 

The solution is a sum of free and forced eccentricities: 



(17) 



const X e 



A 



where 



i-gj 

7 



(18) 



(19) 



The forced solution diverges at secular resonance (7 = 
gj). But this divergence is not physical. It is a conse- 
quence of dropping nonlinear terms. 

With the nonlinear term included, the test particle's 
Hamiltonian become^ 



H{z) = 7 



\zr 
4 



(eje 



igjt . 



+ C.C.) 



with equation of motion 

1 dz 
ij dt 



zil 



) - eje 



igjt 



(20) 



(21) 



The nonlinear term reduces the free frequency from 7 to 



:7(1- 



Y2). We call g the nonlinear free frequency. 



It can also be introduced by first re-writing the Hamil- 
tonian in terms of the real canonical variables Pe and zu 
and defining 



dw I . . . 

-U = .(l-,e/2). 



(22) 



Because ^ is a function of eccentricity, secular resonance 
can occur if the eccentricity is chosen so that g ^ gj. 
This is a nonlinear secular resonance. 

Solutions of equation (21) are shown in Figure [l] for 
two cases, j < gj and j > gj. When j < gj (left panel), 
the free frequency at small \z\ is less than the forcing 
frequency, ^li^i^o ^ 7 < 9 J- Increasing \z\ decreases g, 
moving it further away from resonance. Hence the par- 
ticle's evolution is similar to the linear case, but with a 
smaller precession rate. When j > gj (right panel), the 
precession rate at small values of \z\ exceeds gj. Increas- 
ing \z\ causes the precession rate g to decrease until it 
nearly matches gj. A resonant island appears, within 
which the angle w — gjt librates. Increasing \z\ even fur- 
ther forces the free frequency to be less than ^j, taking 
the particle out of secular resonance. 

To calculate the width of the nonlinear secular reso- 
nance, we first write Hamiltonian (20) in terms of the 
real canonical variables p^ and w (eq. [s]), and then 
convert the Hamiltonian to a time-independent one by 
making the canonical transformation to (pe,'^-), where 



tZ7_ = tu - gjt , 



(23) 



^ We drop the C3 term (see Table [T]) because it merely alters 
the coefficient of the ci term by a smaTTamount. 
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Fig. 1. — Trajectories of z for a test particle in the Presence 
of a Coplanar Jupiter: T he curves are solutions of Hamiltonian 
( |20| , i.e. of equation ( |21[ ). Equivalently, they are level curves of 
Hamiltonian (l24|. The y-axis is approximately e^, and the a:-axis 
is the phase, modulo 27r. Nonlinear secular resonance occurs for 
the case gj < 7. 



yielding the new Hamiltonian 



7 4 



2e 



jWPe COS tU_ 



(24) 



This form for a Hamiltonian has been called the "sec- 
ond fundamental model for resonance," the first being 
the pendulum. Its properties have been extensively cat- 
alogued because it also applies to first-order mean mo- 
tion resonances, for which the variables have a different 



interpretat ion (Henrard & Lemaitre 1983 [Murray fc Der- 
[mott|2QQQ ). The numerical integrations shown in l^lgure 
[T] trace level curves of this Hamiltonian. For parameters 
such that the resonant island both exists and is suffi- 
ciently far from \z\ =0 (as in the right panel of Fig. 
IT]) , one may approximate the multiplying the cosine 
Term as constant, in which case the Hamiltonian is that 
of the pendulum. The center of the resonant island is 
where {d/dt)w- = 0, i.e. at 

Pe. = 2A (25) 

in this approximation; this is equivalent to g = gj. The 
width of the island may be found by first completing 
the square in the "kinetic" part of the Hamiltonian, i.e. 
setting —p1/4-\-A'Pe = — |(pe — Pe*)^+const. Since H_ 
is constant, the half-width of the island is given by 



Spe 



:4ey'(2A)V^ 



(26) 



3.2. Jupiter and Saturn 

We now add in the effect of a second coplanar planet, 
Saturn, and assume that Jupiter's and Saturn's evolution 
is described by their linear Laplace-Lagrange solution. 
In that solution, Jupiter and Saturn participate in two 
normal modes, which we call the Jupiter-dominated and 
Saturn-dominated modes. We denote the eigenfrequen- 
cies of these two modes gj and gs- Jupiter's (complex) 
eccentricity is a sum of two terms, one for each mode. 



and may be written as ej^je 



igjt 



ej,se 



igst 



Similarly, 



Saturn's eccentricity has one term cx e*^^^ and another 
(X e*^-^^. When all four terms are included, one arrives 



w - gjt 



Fig. 2. 
eq. M) 



- Surfaces of Section for the Coplanar Case {H given by- 
Section taken at times when e^^^J~^s)i = i- the y-axis 
is approximately e"^; the a:-axis is the phase, modulo 27r. In the left 
panel, the chosen parameters yield non-overlapping separatrices, 
and very little chaos is seen. In the right panel gj is slightly larger, 
yielding overlapping separatrices and a sea of chaos. 

at the following form for the test particle's Hamiltonian 
(see eq. 20 



H{z) = 7 



(eje 



igjt . 



- c.c.) 
(27) 

where 7 is now re-interpreted to represent the test parti- 
cle's precession rate due to both Jupiter and Saturn. The 
ej term is due to the Jupiter-dominated mode, with cj 
now a weighted sum of both Saturn's and Jupiter's eccen- 
tricity within that mode; similarly, is for the Saturn- 
dominated mode. For the purposes of this section, rather 
than solving the linear Laplace-Lagrange equations, we 
consider the parameters 7, ej, and es to be adjustable 
constants. 



Hamiltonian (27) would also describe the test parti- 
cle's evolution if Jupiter's and Saturn's eccentricities and 
precession rates were enforced by hand to be 67,65,^7, 
and gs^ with ej = \ajej and es = ^(^s^s- This re- 
interpretation of Hamiltonian (27), while unphysical as 



far as Jupiter and Saturn are concerned, can be helpful 
when considering the test particle's evolution under their 
influence. 

Figure [2] shows results of numerical integrations of 
Hamilton's equation, plotted as surfaces of section. 
Changing to real canonical variables (eq. [s]), we see 
that Hamiltonian (27) has two cosine terms. The ej 
cosine term acting alone would yield a resonant island 
as long as 7 > ^j. Similarly, the term would yield 
an island if 7 > gs- The location and width of the is- 
lands are quantified in Figure |T] When acting together, 
there may be two resonant islands. The left panel of 
Figure [2] shows a case when the parameters have been 
chosen to yield two non-overlapping islands. The result 
is mostly regular motion. The right panel shows what 
happens when gj is increased, so that Jupiter's A (eq. 
19 ) is reduced sufficiently that the islands overlap: the 
overlapping islands break up into a sea of chaos. This 
the we ll-known Chirik ov resonance-overlap criterion for 
chaos d Chirikov 1979). From the widths and locations 
of the resonances as displayed in Figures [l] and |2j one 
deduces that the criterion for chaos is 2\gj — gs\h ^ 
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46V^(2A)V\+46y ^(2A.c,)V4^ where As ^ {l - gs)h 



1/2 



( |SidhchQvsk^[T99Q] ). 

it is instructive to consider the form of the surface of 
section shown in Figure |2] in more detail. This reason- 
ing will also be helpful when we include a second degree 
of freedom (i.e. inclination) below. When the motion 
of z is regular, it can be written as a Fourier sum of 
terms with frequencies equal to the three fundamental 
frequencies of the problem (i.e. g^gj^gs^ where g is the 
nonlinear free frequency) , as well as integer combinations 
of these frequencies. But since only relative frequencies 
are physically meaningful, there are really only two fun- 
damental frequencies, which may be chosen to be ^ — gj 
and gs — gj, i.e. relative to gj. In other words, if z is 



regular, then z = ze 



-igjt 



is a doubly periodic function 



with periods 2ti /{g — gj) and 27r / {gs — gj)- In Figure [2j 
we choose to plot the amplitude versus phase of z when- 
ever the second period completed an integer number of 
cycles. As long as z is regular, the value of z at those 
times is a singly periodic function, and hence appears in 
the plot as a connected curve. By contrast, when z is 
chaotic, it appears as scattered points. 

4. AN INCLINED AND ECCENTRIC JUPITER 
4.1. Equations of Motion 

In this section, we consider the evolution of a test par- 
ticle that comes under the influence of a single planet 
("Jupiter") that has fixed values of eccentricity ej, in- 
clination zj, apsidal precession rate gj = zuj^ and nodal 
precession rate sj = Ctj. 

This is a model for the case when a particle comes 
under the infiuence of two planetary Laplace-Lagrange 
modes, one eccentric and one inclined. (See the discus- 
sion of Hamiltonian [27].) As we show below, in the 
real Solar system the main modes affecting Mercury are 
the Jupiter-dominated eccentricity mode and the Venus- 
dominated inclination mode, with the Venus-dominated 
eccentricity mode also playing a role. Therefore for appli- 
cation to Mercury, ej and gj refer to the amplitude and 
frequency of the Jupiter-dominated eccentricity mode, 
while ij and sj refer to those of the Veni^^-dominated 
inclination mode. Nonetheless, for the purposes of the 
present section it is simplest to assume that Jupiter is 
the only planet, and that both of its precession rates are 
enforced by divine intervention. 

We evolve the secular equations for the test particle to 
leading nonlinear order. Because the test particle now 
has two coupled degrees of freedom, its evolution is more 
complicated than before, and there are man y more t erms 
to include in its Hamiltonian. At first (§^ 41|[X4 ), we 
include only the following terms from Table [Tj 



^-H{zX) = \z\'-\C\'-^-^ 



\C\' 



2\z\'\C\' 



-{eje'^'h* - ije'''\* + c.c.) , (28) 

where ej = |aej. Note that the effect of Jupiter's ec- 
centricity is diluted by a factor ~ a <C 1, whereas the 
effe ct of its inclination is undiluted by any such factor. 
In ^4.5 we add in all the remaining terms from Table 
|T] and show that these additional terms have little ef- 
fect in the parameter regime of interest. The terms in 
Hamiltonian ( 28 ) that are second order in eccentricity or 



inclination (i.e., the first two terms, and the bracketed 
terms on the second line) are responsible for linear evo- 
lution. We keep only three nonlinear terms, oc |2:|^, 
and l^plCp. As we show in this subsection, these are re- 
sponsible for nonlinear frequency shifts. And as we show 
in subsequent subsections, nonlinear frequency shifts are 
crucial for resonance overlap and chaos. Even though the 
frequency shifts might be small (second order in eccen- 
tricity and inclination), they can still be sufficient to shift 
the frequency into and out of secular resonance. One of 
the terms we drop is the Kozai resonance, i.e., the term 



C26 X 

feet on 



+ C.C.) in Table jT[ That term has little ef- 
1 Mercury's evolution. Because its phase is rapidly 
varying, and hence the term near ly av erages to zero for 
parameters similar to Mercury's ( ^4.5[ ). By contrast, the 
frequency-changing terms can never average to zeroj^ 
The equations of motion are 



1 dz 
ij dt 
1 dC 
ij dt 



-z{l- 



\zl_ 
2 

2 



2|CH - eje'sjt 
-2\z\^) + ije'''* . 



(29) 
(30) 



Expressed in terms of the real canonical variables (eq. 
[12]), Hamiltonian ( [28| ) has two cosine terms (i.e., pri- 
mary resonances), wmch become important when their 
arguments vary slowly. We first focus on the term 
oc cos(tz7 — gjt). Defining the free nonlinear apsidal fre- 
quency 



dw I 

dt ^^J=^J= 



(31) 



one would expect this resonance to be important near 
where 

g{Pe^Pi) = gj^ (32) 



as may also be inferred by an inspection of equation (29). 
To nonlinear order, ^ is a function of eccentricity and 
inclination, because of the nonlinear frequency-changing 
terms included in Hamiltonian (28). Hence by varying 
e and z, one can alter the free precession frequency and 
bring it into secular resonance. In the Pe-Pi plane, the 
resonance traces out a one dimensional curve — or, in fact, 
a straight line to leading nonlinear order. 

The second cosine term in the above Hamiltonian, 
cos(l] — sjt), behaves similarly. Defining 



dn^ 



■1 



1 



^Pi - ^Pe 



(33) 



one would expect it to be important near where 

s{pe,Pi) = sj. (34) 



We shall make these considerations more precise in §4.3[ 
where we also work out the resonant widths, and show 
that in addition to the primary resonances are a multi- 
tude of secondary resonances. 

^ We do drop some frequency-changing terms, specifically ones 
that are given by const. x|2;p, and const. x|("p, where the constant 
is 0(e^,i^). Even though these terms do not average to zero, they 
merely shift the linear frequencies, and hence do not change the 
behavior qualitatively. In the absence of these terms, the linear 
apsidal and nodal frequencies are equal and oppo site ; we rectify 
this shortcoming in our ^c-model Hamiltonian (eq. [53]). 
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Fig. 3. — Map of the Mean Momenta (MMM) with low excitation: 
Each point is the result of a single integration of Hamiltonian ( |28| , 
with parameters ej = ij = .0003, gj = 0.727, and sj = — 1.2(37. 
The y-axis is the time-averaged jCj ~ ^ind the x-axis is the time- 
averaged ~ e^. The initial conditions are on a uniform grid in 
jCp and jzp. Resonant bands are clearly visible, as is the zone of 
chaos where the bands overlap. Because of the averaging, chaotic 
orbits that explore the full extent of the overlap zone give rise to 
points at the center of the zone. Chaotic orbits that only partially 
explore the overlap zone (as is true of the real Mercury today; see 
below) give rise to chaotic points at the edges of the overlap zone. 
The orbits are very regular at small pe, Pi] for comparison, the real 
Mercury currently has {pe) ~ .05 and (pi) ~ .02. Its motion would 
be regular under this weak forcing. 



4.2. Simulations: Maps of the Mean Momenta 



We run suites of simulations of Hamiltonian (28), i.e., 



of equations (29)- (30). There are only five parameters, 
7, ^j, sj, ej, and zj. The linear frequency 7 sets the 
overall timescale, and can be scaled out. We choose 
gj/j = 0.72 and sj/7 = —1.26, since these are close 
to the true values in the Solar system for the Jupiter- 
dominated eccentricity mode (relative to Mercury's free 
precession frequency), and for the Venus-dominated in- 
clination mode (see ^ . We also set the excitation ampli- 
tudes equal to each other, ej = zj, and present sequences 
of simulations with various amplitudes. The true value 
in the Solar system for the corresponding modes is, very 
roughly, ej ^ ij ^ 0.01 (see ^for more precise values). 

The dynamics is that of two nonlinearly coupled har- 
monic oscillators, each of which is also nonlinear and 
is forced periodically. We have attempted many differ- 
ent methods for visualizing the integration results, such 
as using catalogs of surfaces of section or Fourier trans- 
forms. However, most methods were too complicated, 
and obscured the underlying simplicity of the dynamics, 
i.e. that it is the overlapping of resonances that drive 
chaos. In the end, we invented a new method, the Map 
of the Mean Momenta (MMM). This method has many 
advantages over the usual surfaces of section. It is some- 
what similar to the frequency map analysis of iLaskar, 



( |1990D (see below). 

iTgure|3]maps the results from around 50, 000 numeri- 
cal integrations at very small excitation, ej = ij = .0003, 
using the MMM. Each point in the plot is the time- 
averaged value of |2;p ~ and |CP ^ i'^ from a single 
simulation, averaged over a time span of 3 x 10^/7. Be- 
fore t aking the tim e average, we filter with a Hanning 
filter ( Laskar|[l993 ), which leads to a faster convergence 
of the averages (when they do converge). The simula- 
tions were initialized with values of |zp and |CP that 
were equally spaced on a grid, with the spacing in \z\'^ 
twice that in the initial phases were tu = 7r/2 and 
n = -7r/2. 

Three kinds of motion are readily apparent in the 
MMM: (i) regular and non-resonant, (ii) regular and res- 
onant, and (iii) chaotic. Most of the figure is covered 
with a regular grid of points that nearly traces the initial 
conditions. Here, the values of z and ( remain regu- 
lar and non-resonant throughout the simulation. A few 
resonant bands also appear in the figure, where the mo- 
tion is also regular. Note that if the initial conditions 
span a resonant island, and if the motion remains regular, 
then the time-averaged momenta (|zp and |CP) exhibit 
a sharp discontinuity. Inside the island, they average to 
their values near the island center, whereas outside the 
island's separatrix they average to a value offset from the 
center by a finite amount, of order the resonance width 
(or, more accurately, around 1/4 of the resonance full- 
width). This can be seen clearly in Figure Isj where the 
regular points at the center of the resonant oands repre- 
sent librating particles. Also apparent in the figure are 
the chaotic trajectories. These show up as the cluster of 
irregular points near where resonant bands intersect. We 
have checked the Lyapunov exponent, as well as surfaces 
of section (see below), to verify that points that appear 
on the figure to be chaotic are truly chaotic. 

Figures [4][5] show the MMM for simulations with the 
same parameters as in Figure |3j but with ej and ij in- 
creased first to 0.003, and then to 0.01. With increas- 
ing forcing, the locations of the resonant bands do not 
change, but they get wider, and higher order resonances 
become visible. As a result, the zone of chaos where 
the bands overlap expands. Surprisingly, even with the 
seemingly modest forcing of ej = ij = 0.01 — values that 
are comparable to those in the real Solar System (see 
below) — the zone of chaos approaches very low values of 
e and and close to the values for the real Mercury. 

Our method for displaying results, the MMM, is similar 
in ph ilosophy to frequency map analysis (FMA; Laskar 
1990 ). But whereas in FMA one plots the frequencies 
of the co-ordinates, here we plot the averages of the mo- 
menta. We have also performed the FMA (not shown); 
and when we convert from frequencies to momenta via 
the inverse of equations (31 ) and (33), the resulting maps 
are almost identical to theMMM. Tor the purposes of the 
present paper, we prefer the MMM because its axes are 
approximately (e^) and (z^), which are simpler to inter- 
pret than the precession frequencies. 

4.3. Theory: Resonance Locations and Widths, and 
Zone of Chaos 

To develop understanding of the behavior seen in the 
MMM's, we first re- write Hamiltonian (28) in terms of 
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Fig. 4. — MMM with Medium Excitation: Similar to Figure 
but with ej and ij increased by a factor of 10. The resonant bani 
are larger, as is the chaotic zone, which has encroached much closer 
to where the real Mercury lies, ~ (.05, .02). The dashed magenta 
square is for comparison with the axes of Figure [s] 




Fig. 5. — MMM with High Excitation: Similar to Figures[3p] but 
with ej = ij = 0.01. Note the expanded scale. Many high order 
resonances are visible. The zone of chaos approaches the origin, 
even though ej and ij are <C 1. The points labelled Mercury are 
the result of an N-body simulation of the full Solar system; each 
point is Mercury's pe and pi averaged over a timespan of 100 Myr, 
for the first 600 Myr of the simulation shown in Fig. |16| Mercury's 
true orbit lies near the boundary between regular motion and chaos 
in the MMM of the simplified model. Parameters used for this 
MMM are within ~ 20% of the the true Solar system values. The 
true Solar system is more chaotic due to other forcings. 



the real canonical variables (eqs. [8]-[9]), and then trans- 
form from (pe, ^) to (pe, ^- = ^ — gjt) and from (p^, 
to (j9f , = Q — sjt), which transforms the Hamiltonian 
to a time-independent one: 



7 



■Pi 



Pe - Pi 



—2pePi — 2ejy^costz7_ + 2zjy^cosl]_ , (35) 

where A is the linear apsidal frequency mismatch (eq. 
[19]), and 



A. = 



-7 - sj 
7 



(36) 



is the mismatch for the nodal frequencies. In the absence 
of the coupling term (ex PePi) both degrees of freedom 
would evolve independently accord ing to the equations 
of the second fundamental model (^3.1). 

In the following, we determine the location and width 
of the two primary resonances: the eccentricity resonance 
([1,0]) and the inclination resonance ([0,1]). To do so, we 
ignore inclination forcing (setting ij = 0) when studying 
the eccentricity resonance, and vice versa. The system 
is trivially integrable if either zj = Oorej = 0. In the 
former case, Pi is constant because the Hamiltonian does 
not depend on Therefore H- is equivalent to the 

coplanar Hamiltonian (eq. [24]), but with A ^ A — 2pi. 
This implies that the center or the eccentricity resonance 
is located at 



Pe 



2(A-2p,) ^ [1,0] 



(37) 



(eq. [25] ) and the island has half- width at fixed Pi given 
by 

^Pe = 4(ejv^)i/2 ^[1,0] (38) 

(eq. [26]). The interior of this island is plotted in Figure 
[6] as aolue band, with parameter values as chosen for 
the simulations of Figure [3] Comparing the two figures 
shows that the above analytic expressions agree well with 
the result of the numerical integrations. Note that we 
plot the half- width in Figure [6] rather than the full- width 
because the average momentum of an orbit that lies just 
outside of a separatrix is approximately half-way to the 
edge of the resonance. Figure [6] also shows the inclination 
resonance as a red band. Reasoning as before, its center 
and half-width (at fixed Pe) are 




[0,1] 
[0,1] 



(39) 
(40) 



In addition to these primary resonances are an infinite 
number of secondary ones. Far from resonances, we may 
ignore the cosine forcing terms since they tend to average 
to zero. The test particle's free precession frequencies 
relative to Jupiter are then 

g-^g-gj^ ^ L.^^.^o = 7 ( A - - 2^,141) 



dt ^^J=^J- 



1 



Resonances are important near where 

mg- + ns- =0 [m, n] 



(43) 
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Fig. 6. — Four strongest resonances and their widths and re- 
gion of overlap (analytic calculation): Parameters are as in Figure 
p] The center of each [m, ril resonance is the line determined by 
mg- -\-ns- = (eqs. ^^-p2]). The arrows show the half-widths, 
with the orientation aligned with the direction of motion in the 
corresponding resonances. The central grey shaded region shows 
the region of overlapping separatrices, where chaotic motion is ex- 
pected. Comparing with Figure [3] shows agreement between theory 
and simulation. 

for integer pair [m, n]. Therefore the center of each [m, n] 
resonance traces out a Une in the Pe-Pi plane. For in- 
stance, the centers of the [1,0] and [0,1] resonances are 
as worked out above, and the center of the [1,-1] reso- 
nance is the hue 

Pe* - lp^. = ^ (A, - A) ^ [1, -1] , (44) 

which passes close to the origin for our choices of A and 
Ag. In general, the slope of a resonance line in the Pe-Pi 
plane is (4n + m)/(n — 4m), and all [m, n] resonant lines 
intersect at the point in the Pe-Pi plane where g- = S- = 
0, i.e. at (pe**,Pi**), where 

2 2 

Pe** = ^(A + 4A,), pi^^ = ^(4A - A,) (45) 

Although there are an infinite number of secondary 
resonances, most are very weak, i.e. their widths are 
small. The most prominent resonances in Figure |3] are 
the primary resonances [1,0] and [0,1], whose widths have 
been worked out above. Next most prominent are the 
[1, -1] and [1, 1] resonances, whose widths may be un- 
derstood qualitatively as follows (see Appendix B for a 
quantitative calculation). The linear solution for z is 
a sum of two terms, the free and forced complex ec- 
centricities (eq. [is])- Similarly, to linear order ( is 
a sum of free and forced complex inclinations. There- 
fore to leading nonlinear order, the coupling term in 
the Hamiltonian {pePi = I^PlCP) can be written as a 
sum of terms, one of which has the form z^jyZ^^^Q ~ 
6(1)6 fi(f)i f ex.-p{i{g(f) — gj — S(f) -\- sj)), where the subscript 



/ denotes forced and (j) denotes free. This term has fre- 
quency corresponding to the [1,-1] resonance; hence the 
width of this resonance is ^ y/\e^i^efif\. The [1, 1] 
resonance behaves similarly. 

The quantitative calculation in the appendix shows 
that, in agreement with the abov e est imate, the [1, ±1] 
resonances have half- widths (eq. [B10| ) 

6pe = Spi = 2^ 6(1)1(1)6 fif [1, ±1] , (46) 
after defining the free values as 

and the forced values as 

^ \A-Pe./2-2pi,y ^ |A, +p^*/2-2pe*| ^^^^ 

where the asterisk denotes values at resonant center, and 
we neglect here the small difference between the lower- 
case momenta and the upper-case ones used in the ap- 
pendix. As in linear theory, the forced eccentricity scales 
as the inverse of the frequency detuning (eq. [iS]), al- 
though now it is the nonlinear frequency detunmg g- 
(eq. [41]) that is relevant; similarly, ^/ is inversely pro- 
portional to 5_. 

The above widths for the [1,±1] resonances are shown 
in Figure |6l The [1, -1] resonance leaves Pe + Pi nearly 
constant, which produces trajectories in the Pe-Pi plane 
that have a slope of -1 (see Appendix B; we again neglect 
the small difference between lower- and upper-case mo- 
menta). Therefore for each value of (pe*:Pi*) at resonant 
center (i.e., where g- — S- =0), the upper and lower en- 
velopes of the resonant band are at (pe* ±^Pe,Pi* T^Pe)- 
Comparing with Figure |3] shows excellent agreement be- 
tween theory and simulations. 

The central grey region in Figure [6] is the overlap zone 
of the above four resonances. Its shape is peculiar be- 
cause each resonance induces a trajectory with a par- 
ticular orientation in the Pe-Pi plane, and we take the 
overlapping region to be wherever one resonance can in- 
duce motion into a second resonance. The zone of chaos 
seen in Figure^ is about twice as large as that predicted 
in Figure |6] This difference is not surprising because we 
plot the half- widths in Figure [6) whereas one might ex- 
pect that chaos would begin where the full-widths over- 
lap. Furthermore, we shall show that higher order com- 
binations of the four strongest resonances also play a role 
in the chaos. Nonetheless, the grey zone provides a rea- 
sonable estimate of the zone of chaos. 

Figure [7| plots the widths calculated above when the 
forcing ej and ij take on the values of Figure [4] The 
latter figure is overlayed on Figure [7| , showing excellent 
agreement between theory and simulation in the zones of 
regular motion. And, as before, the zone of chaos from 
the simulations is around twice as large as the region 
where the resonant half- widths overlap. 

The dimensions of the zone of chaos may be estimated 
analytically. These estimates will be used later to infer 
the threshold for chaos in the Solar system. In Figures 
[6][7| the horizontal and vertical spikes of the chaotic zone 
are due to the [1, 0] and [0, 1] resonances; the extent 
of the spikes is simply estimated by the resonance width 
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Fig. 7. — Resonant bands and their region of overlap, over- 
layed on MMM: The colored bands show the theoretically predicted 
widths of the four strongest resonances, as in Figure [6] The grey 
region is the zone of resonance overlap for these four resonances. 
Also shown is the MMM of Figure |4] showing excellent agreement 
between the theoretical and numerical resonant widths, and satis- 
factory agreement for the zone of chaos. 

near zone center, 



4e 



1/2 1/4 



Pe 



^ [1,0] X [0,1] (49) 
^[1,0] X [0,1] , (50) 



for the horizontal and vertical spikes respectively, where 
Vf^^^ and Pi^^ are the co-ordinates of the zone center (eq. 
[45]). The extent of the zone caused by the overlap be- 
tween the [1,-1] separatrix with the [1,0] resonance can 
be estimated similarly as 

^Pe* = ^Mj)'/^(Pe**Pi**)'/' ^ [1,-1] X [1,0] (51) 

where k (20/17) y^3/4 is an order-unity constant rl The 
expression for 6pi^ is the same, as are the extents oue to 
the overlap between either [1,1] or [1,-1] with either [1,0] 
or [0,1], albeit all have different order-unity values for k. 

It might appear surprising that the width of the chaotic 
zone due to the overlap of a primary resonance with a 
secondary one (eq. I5l]) is not much smaller than that 
due to the overlap between two primary resonances, de- 
spite the fact that the width of a primary resonance is 



Fig. 8. — Su rfaces of Section from integrations of Hamiltonian in 
equation ( |28| >: The parameters are gj = 0.727 and sj = — 1.267, 
and other parameters as shown. The three upper left panels show 
surfaces of section with increasing forcing, showing how the [1, -1] 
resonance gets wider, and its separatrix breaks up into a sea of 
chaos. The (blue) chaotic trajectory in the upper right panel has, 
very roughly, parameters comparable to the real Mercury. The 
lower right panel shows a "double section" of this blue chaotic tra- 
jectory. In the double section, it traces out a branch of a hyperbola 
in momentum space. Comparing with Figure [s] shows where this 
trajectory lies relative to the MMM. 

eccentricity and inclination (eq. [48]). Because of this, 
the [1,±1] resonances play an important role in setting 
the extent of the chaotic zone. For example, with the 
forcing frequencies that we have chosen for the MMM's 
(which are comparable to those for Mercury), the [1,-1] 
resonance overlap allows the region of chaos to encroach 
upon the origin (e ~ 0, z ~ 0) at lower values of the 
forcing than would have been expected based solely on 
the overlap of the [1, 0] and [0, 1] resonances. In ^ 
we shall show explicitly that the [1,-1] resonance plays a 
dominant role in driving chaos for the real Mercury. 

One could also proceed to calculate the widths of 
higher order resonances. Extending our reasoning from 
above, the width of an [m, n] resonance in the Pe-Pi plane 
should scale as 



5p ' 



'O(el^lzl^l) 



[m, n\ 



(52) 



38 



40 



first order in eccentricity or inclination (eqs. 
whereas the width of the [1,-1] is second order (eq 
The reason for this is that the [1,-1] resonance is en- 
hanced by the denominator that appears in the forced 



^ Focusing on the region to the lower left of the zone center in 
Figure [6] or m the center of the [1,-1] resonance is displaced from 
zone center m the Pe-Pi plane by the vector —(l,3/5)x (eq. [44] ); 
we take x > 0. The orbit given by the half-width of the [1,-lJ is 
therefore displaced from zone center by the vector — (1, 3/ 5)x + 
(-l,l)(5pe, where 6pe - 2(2?e**25^**)^/^yw/(17x/10) (eq. 
Equating that vector to the displacement of the center of th e |l'0] 
resonance, i.e. to (— l,l/4)a;^ (eq. [37]) yields equation ( |5l| . 



where e is comparable to the typical eccentricity, and 
i is comparable to the typical inclination. However, as 
we have seen for the [1,±1] resonances, near-resonant 
denominators can make the widths significantly larger 
than this naive estimate. 

4.4. Surfaces of Section and Libration of Resonant 
Angles 

Here, instead of the global map (MMM), we investigate 
a few particular trajectories in detail to demonstrate the 
chaotic behavior, using both surfaces of section and res- 
onant angles. 
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Figure |8] shows a number of surfaces of section from 
Hamiltonian (28). Since the transformed Hamiltonian 
(eq. [35]) is time independent with two degrees of free- 
dom, we fohow the usual practice of taking a section 
whenever the phase of one of the degrees of freedom 
(here, ti7_) executes an integer number of cycles. At 
those times, we plot the amplitude versus phase of the 
second degree of freedom {pi vs. rt_). We may also 
understand this form for the surface of section as fol- 
lows. Hamiltonian (28) has four fundamental frequen- 
cies, ^j, and s wnere g and s are the nonlinear free 
frequencies of z and respectively. Only relative fre- 
quencies are physically meaningful, and there are three 
of these, which we may choose tohe g — gj, s — sj, and 
g — s. But Hamiltonian (28) does not depend on — 1], 

s does not enterj^ There- 
fundamental frequencies, as 
And, as described there. 



and hence the frequency g - 
fore there are two remaining 
for the coplanar case ( ^3.2| ). 
to examine the characteristics of the motion, one may 
take a section whenever the phase corresponding to one 
of the fundamental frequencies (here, the phase vu — gjt 
which corresponds to g — gj) executes an integer number 
of cycles. 

The three upper left panels of Figure |8] show surfaces of 
section with values of gj and 5 j as before, and with vari- 
ous values of ej and ij. All the surfaces of section shown 
have the same (constant) value of energy, H_ = 0.01 547. 
To map out phase space would require many different 
energy values, but for that purpose the MMM is more 
useful. The lower right panel in Figure |8] shows a "dou- 
ble section" of the (blue) chaotic trajectory in the upper 
right panel, i.e. it shows the two momenta wherever both 
vu- and have executed a half-integer number of cy- 
cles. At these times, the cosine terms in the transformed 
Hamiltonian vanish, and all trajectories with a fixed en- 
ergy H_ fall along the same branch of a hyperbola. Of 
course double sections from all the trajectories shown in 
Figure [S] would lie along the same branch of the hyper- 
bola because they all have the same value of H_. 

The surface of section in the upper left panel has a rel- 
atively low €j and ij, and the motion is mostly regular 
for the value of energy chosen. This can also be seen in 
the MMM (Fig. [4| near the relevant hyperbola branch. 
The [1, -1] resonance is clearly evident in the top left 
panel of Figure |8] Its half- width is 5pi = 0.0055, as com- 
pared to the prediction of 0.006 from equation (46). The 
upper right panel of Figure fs] shows the case with higher 
forcing. With this higher forcing, the [1, -1] resonance 
is wider, and the region near its separatrix has broken 
up into a wide zone of chaos. Some of the higher order 
resonances are visible in this section. In the lower left 
panel, the forcing has been raised further. Even though 
ej and ij are still relatively small compared to unity, the 
zone of chaos is vast. 

The blue chaotic trajectory in the upper-right panel 
of Figure [S] behaves qualitatively like the real Mercury, 
and the parameters are also similar (^. Therefore we 
investigate it in more detail. From its surface of sec- 
tion, we see that the separatrix of the [1, -1] resonance 

The full fourth order Hamiltonian does depend on vj — Q 
because of the Kozai term; see Table ^ Therefore one cannot take 
surfaces of section of the fu ll fo urth order Hamiltonian, but one 
can stih plot its MMM (Fig. [Tol. 
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Fig. 9. — Chain of librating angles: Each panel shows resonant 
angles mm- + nfl- (modulo IOtt) with various values of [m,n], 
for the blue chaotic trajectory of Figure [s] The green shaded 
zones show librating angles. Different resonant angle combina- 
tions librate in turn. The first grey strip is when [6,-5] librates 
(not shown). The second grey shaded strip shows a time when the 
[1,-1] and [4,-5] alternately librate in rapid succession, and no other 
angles are clearly librating. 

is largely responsible for driving the chaos for this or- 
bit, together with overlapping higher order resonances. 
This orbit remains bounded by the [3,-2] and [2,-3] res- 
onances, and hence can never come under the direct in- 
fluence of the primary resonances ([1, 0] and [0,1]). The 
bound on the chaotic zone is a consequence of using a 
truncated Hamiltonian that can be written in a time- 
independent form with two degrees of freedom. For the 
full Hamiltonian, one might expect that diffusion could 
act on long timescales (Arnold diffusion), ultimately al- 
lowing the trajectory to cross into other regions of phase 
space. 

Figure [9] shows explicitly that the chaos is due to the 
overlapping of high order resonances. The resonant an- 
gles mvu- + nft- are plotted for various values of [m, n]. 
Different resonant angles librate in turn, showing that 
this orbit first comes under the influence of the [1,-1] 5 
then the [5,-4], then the [1,-1] 7 etc. 

4.5. Full Fourth Order Hamiltonian 

Thus far we have focused on the truncated fourth order 
Hamiltonian (eq. 



28 ). Figure 10 shows the MMM of the 



full fourth order Hamiltonian, expanded to leading order 
in a (i.e., including all terms in Table fl] in Appendix A). 
From the similarity of Figure [lO] to trie truncated inte- 
grations of Figure [5) we conclude that the terms dropped 
in the truncated Hamiltonian have little effect on the dy- 
namics, particularly in the region of small e and i (lower 
left corner of the MMM). 

The dropped terms have little effect because the only 
new fundamental relative frequ ency they introduce is ^ — 
s (see first paragraph of ^.4). This "Kozai frequency" 
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Fig. 10. — MMM with full fourth order Hamiltonian: Similar 
to Figure [s] but integrations have been performed with the full 
Hamiltonian (Ta blefTTin Appendix A), rather than the truncated 
Hamiltonian ( |28| . From the fact that the two figures are broadly 
simi lar, one can infer that the terms dropped from Hamiltonian 
(128} are of small importance in the regime of interest. 



differs significantly from zero in the domain of Figure 10 
and hence it can only combine with the other two relative 
frequencies {g — gj and s — sj) to give resonances at high 
order. To be quantitative, the resonant line of the Kozai 
frequency is at ^ — 5 ~ 7(2 + (3/2)pe — (5/2)pi) ~ 0, 
i.e., it traces the line ~ (4 + 3pe)/5. Hence the Kozai 
resonance is at much larger pi than shown in Figure [lO] 

Aside from the Kozai term (c26 in Table [l]) , all other 
terms dropped from the truncated Hamiltonian depend 
on Jupiter's eccentricity or inclination. These do not 
introduce new forcing frequencies because Jupiter's fre- 
quencies already appear in the test particle's orbit at 
linear order — in its forced e and i. While the dropped 
terms do change the amplitudes of the forcing terms, the 
change is small as long as Jupiter's e and i is smaller 
than the e and i it linearly forces in the test particle, as 
is true of Figure [lO] 

We suspect that the terms dropped from Hamiltonian 
(28) are quite often of secondary importance. This is 



largely true for the real Mercury (Q. And we suspect 
that it is true more generally because if secular interac- 
tions between two planets are strong, then the forced e's 
and z's will typically (though not always) be larger than 
the forcing ones. Nonetheless, the dropped terms can 
be important in certain circumstances; for example, the 
Kozai term will play a role if a planet has a high incli- 
nation, and MMR's will be important for planets whose 
orbital periods are near integer ratios. 

4.6. Fourier Transforms 

In ^we shall make the connection to the real Mercury. 
For that purpose, it will prove instructive to examine 
trajectories in Fourier space. 

For a more exact comparison to Mercury, we consider 
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Fig. 11. — Fourier transforms of th e test particle's z and C for 
the K-model with k = 0.75 (eq. [53| ): With this relatively small 
value of K, the trajectory is quasiperiodic, as indicated by narrow 
spikes in the Fourier transform. The free and forced z are peaks at 
qm and gj in the top panel. The other peaks are due to nonlinear 
couplings, and are at frequencies qm + rn{gM — gj) + n{sM — sv), 
labelled [m,n]. The horizontal red arrows in both panels denote 

frequencies spaced by qm — 9 J- The bottom panel shows that |C| 
is similar, with free and forced C 3,t frequencies sm and sy, and 
nonlinearly generated peaks at sm + ^(9m ~ 9j) + n(sM ~ sy)- 



here the Hamiltonian 

Ih{zX) = \z?-i\C? 

7 



\c\' 



2\zm 



{eje'^^'z"" -ive^'^'C ^c.c.) , (53) 



which differs from Hamiltonian (28) by the inclusion of 



a constant 7 to allow the linear apsidal and nodal pre- 
cession rates to differ from each other (see footnote [8|. 
Note that we also change notation so that iy and sy are 
the amplitude and precession rate of the Venus mode. 
We focus on a one dimensional family of systems pa- 
rameterized by which scales all eccentricities and in- 
clinations. More precisely, in this "A>:-model" we choose 
the parameters ej = iy = 0.008/^; and initial conditions 
\z\ = 0.16/^, Id = 0.07/^, m = Q = 7r/2. The remaining 
parameters are 7 = 0.9, gj = O.727, sy = — I.I47. With 
these parameters, the center of the [1,0] resonance is at 
Pe* = 2(0.28 — 2pi), as before; and the center of the [0,1] 
resonance is at = —2(0.24 — 2pe), whereas before the 
constant was 0.26 rather than 0.24. This difference is 
of little consequence. For displaying the results of the 
integration, we shall choose 7 = 5.87'Yyr. Our ratio- 
nale for choosing these particular numerical values will 
be explained in ^ 

Figure pT] shows the Fourier transforms of z and ( for 
the /t:- model trajectory that has = 0.75. We normalize 
the Fourier transform of z as 



z{uj) 



z{t)e-'"'dt 



(54) 
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and similarly for where T is the duration of the Fourier 
transform, which we choose in the present subsection to 
be T = 4400/7. With this normalization, if z has con- 
stant amplitude and frequency, i.e. if z{t) = /coe*^°^ then 
\z\ = \ko\ at frequencies close to ujq. 

The trajectory used for Figure [TT] is quasiperiodic — the 
peaks in the Fourier transform are simply spikes, whose 
widths become narrower for larger T. We call the two 
largest peaks in the top panel the forced and free z. The 
forced z is at frequency gj = 0.727 = 4.23'Yyr. The 
free z is at frequency gM = 5.78'Yyr. Because of nonlin- 
earities, the free frequency (^m) differs from the linear 
free frequency (7) by a small but non- negligible amount. 
Similarly, in the bottom panel the largest two peaks are 
the forced ( at frequency sy = — I.I47, and the free ( at 
frequency sm = — 5.45'Yyr, which differs from the linear 
free frequency —77. 

In addition to the free and forced z and there are 
a multitude of peaks in Figure 11 that are generated by 
nonlinear couplings. The peaks in z all fall at frequen- 
cies gM + rn{gM - Qj) + n{sM - sy) for integers m, n. 
Roughly speaking, the peak amplitudes become smaller 
for larger values of |m| and |n|. These amplitudes can 
be calculated perturbatively, as is sketched in the fol- 
lowing. As before, we define the free and forced com- 
ponents as z/, C/)^ which have phases that rotate 
with frequencies (^m, -^m, ^y), respectively. To lead- 
ing nonlinear order, the nonlinear terms in the equation 
for dz/dt are proportional to z\z\^ ^A?4> 



/c= 1 



/c= 1.3 /c= 1.55 
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). These generate 

gM),gM ± 



and^lCP = (^0 + ^/)|C0 + C/P (eq. 
six new frequencies in z: (2^m — djJTi'^dJ 
{sm — sv),gj ± {sm — sv)- Each frequency-generating 
term acts as a linear forcing on z. Together with the 
free and forced these account for eight of the peaks 
marked in the top panel of Figure [TTj specifically, they 
account for the two highest peaks in each of the four 
left-most triplets. The other peaks are accounted for by 
higher order nonlinear terms. One of these other peaks — 
the one labelled [1,-1] — is quite large, even though one 
might naively have expected that it would be smaller be- 
cause it enters at a higher nonlinear order. The reason 
for this is that its forcing frequency differs from gM by 
a = {gM — gj) — {sm — sv) which is quite small. Hence 
this near resonance amplifies the peak by gM/cr ~ 20. 
We note parenthetically t hat the width of the [1,-1] I'es- 
onance (as described in ^4.3) is directly related to the 
amplitudes of the three peaks at gM and gM ± cr- The 
Fourier transform of ( behaves similarly to that of with 
the frequency peaks at sm + '^{sm — gj) + ^(^m — sy)- 
Figure [12] shows the Fourier transforms for the />:-model 
at higher values of The left-most panels show the 
case n = 1. The motion is largely quasiperiodic, but the 
nonlinearly generated peaks have increased significantly 
relative to the n = 0.75 case. At n = 1.3 the motion is 
chaotic, and at z^: = 1.55 it is highly chaotic. 

5. MERCURY 

We integrate the eight Solar system planets with the 
SWIFT symplectic integrator (Levison fc D uncan 1994), 
supplemented with a routine lor Me rcury's relativistic 
precession (see 
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We initialize the planets with their current orbits and 
use their actual masses. The integration timestep is 8 
days. 




for code details). ^ seen in Figure 13 



Fig. 12. — Fourier transforms of results from three ^c-model sim- 
ulations, with K = 1,1.3,1.55. At = 1 (left-most panels), the 
motion is still largely quasiperiodic. The amplitudes of the non- 
linearly generated peaks have risen significantly relative to Fig. |11| 
(k, = 0.75), even though the forced and free z and have changed 
by a modest amount. At = 1.3, there is weak chaos — the peaks 
have widened, and neighboring peaks overlap. At = 1.55, the 
trajectory is highly chaotic. 

One might suspect that Mercury's orbital evolution is 
more complicated than our toy model for a variety of 
reasons: its e and i are not too small, and hence the 
fourth order expansion is approximate; it is not massless, 
and hence backreacts onto the other planets (especially 
Venus); there are seven other planets that do not have 
constant orbital elements and frequencies but participate 
in the overall chaos of the Solar system; and Mercury 
can be affected by resonant terms. Despite these com- 
plications, we show that the chaotic behavior of Mercury 
is qualitatively similar to the Hamiltonian model. This 
is perhaps not too surprising, since nonlinear dynam- 
ics are largely driven by resonances and their overlap. 
Hence as long as a model roughly captures the locations 
and widths of the principal resonances, it should produce 
qualitatively correct behavior. 

5.1. Fourier Transforms 

To compare with the /^-model, we first consider cleaner 
cases by pre-multiplying the current eccentricities and in- 
clinations of all planets by the reduction factor A^nbody 
Figure [13] shows Mercury's Fourier transforms in a 
^nbody = 0-75 integration lasting T = ISOMyr. Compar- 
ing this with the model at = 0.75 (Fig. [TT]) shows 
broad agreement. In truth, the parameters lor the k- 
model were chosen to match the free and forced z and 

the frequencies and heights of 



I.e. 



the four peaks marked gj^gM-,sv-,SM- Since there were 
eight quantities to match, we could do this by adjusting 
eight parameters in the /^-model: 7, 7, ^j, 5y , ej, , as 
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Fig. 13. — Fourier transform of Mercury's z and (roughly, its 
complex eccentricity and inclination) in a SWIFT N-body sim- 
ulation, ^Cnbody = 0.75: the initial e's and i's of all planets 
were pre- multiplied by the factor ^Cnbody The peaks here are 
broadly similar to those of the ^c-model. As in Fig. |11| they are 
marked by vertical lines that are displaced from qm and sm by 
^(fi'M — 9j) -\-tt^(sm —sv)- The agreement between the two figures 
shows that the ^c-model captures much of the physics of the real 
Mercury. Nonetheless, there are a number of differences. See text. 



well as the initial values of \z\ and \(\. Therefore it is 
not significant that the forced and free peaks in the two 
figures agree. What is significant is that the other peaks 
that are generated by nonlinear couplings of the forced 
and free peaks also largely agree. This indicates that the 
/^-model captures much of the nonlinearity as seen in the 
real Mercury. 

There are, however, at least three differences of note. 
First, the peaks in Figure 13 are broader than those in 
Figure [TT] This is because Figure 13 suffers from weak 
chaos. But it is remarkable how sharp the largest peaks 
are: even though the e's and z's of the Solar system have 
only been reduced by 25%, the resulting chaos is surpris- 
ingly weak. Note that the integration intervals in the 
two figures are the same, T = 4400/7 = 150Myr, and 
hence the finite width of the peaks in Figure [13] is not 
due to the finite T. A second difference between the two 
figures is that the zm peak at frequency + {qm — 9j) 
is significantly larger in the /^nbody integration. That 
peak is so large because it is overlapped by a peak forced 
by Venus's eccentricity mode, which has precession fre- 
quency gv ^ ^^gu—gj' In other words, for /^nbody = 0.75, 
Mercury is in a secular resonance with a librating angle 
that corresponds to frequency 2gM — gv — gj^ and this 
largely hides the effect of gy in the Fourier transform of 
Figure [13] The third difference between the two figures 




is the peak in C,m that is caused by Uranus's inclination 
mode. But this peak appears to have little dynamical 
consequence for Mercury. 
In Figure [14] the factor multiplying the initial e's and 
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Fig. 14. — Same as Figure [13) but with /^nbody increased from 
0.75 to 0.95. The motion is more chaotic here, and Venus's eccen- 
tricity forcing at frequency gv is distinct. The blue dashed arrows 
denote frequency spacings of gv — gM : a-nd the blue dashed verti- 
cal lines denote peaks due to the gv mode and nonlinear couplings 
generated by that mode. 



i's has been raised to /^nbody = 0.95. The resulting mo- 
tion is more chaotic, as the widths of the peaks are wider 
than before, especially for (m- In addition, the frequen- 
cies have been shifted sufficiently to break Mercury from 
the 2gM — gv — gj resonance, and the effect of Venus's 
eccentricity forcing is distinct. Even though the motion 
is more chaotic, the principal forcing peaks and their 
harmonics are still identifiable. 

Figure [T5] shows the Fourier transform of the real Mer- 
cury (^^nbody = 1)- Even though the initial e's and i's 
have been increased by only 5% relative to Figure [m] 
the motion is significantly more chaotic, and the nonlin- 
early generated peaks are less easily identifiable, espe- 
cially those near Mercury's free frequencies gM and sm- 
Nonetheless, we conclude from the progression of Figures 
[T3][T5l that the /^-model captures much of the physics. 
In particular, two modes — the Jupiter eccentricity mode 
and the Venus inclination mode — are primarily respon- 
sible for driving Mercury's chaos. The most important 
element lacking from the /^-model appears to be the ex- 
tra forcing by the Venus eccentricity mode. The fact 
that the A>:-model becomes chaotic at a higher threshold 
than the real Mercury {k, ~ 1.3 vs. A^nbody ^ 1) is likely 
partly due to that extra forcing (i.e. that extra forcing 
makes the real Mercury more chaotic). An additional 
contributor to the discrepancy between the two critical 
hz^s is that the simple /^-model does not accurately cap- 
ture nonlinear frequency shifts, while the precise values 
of the frequencies are important for where the resonances 
overlap. Despite this, the difference between the critical 
/t:'s is not large, and this lends support to our claimed 
origin for Mercury's chaos. 

We note parenthetically that while we only focus on a 
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Fig. 15. — Same as Figs. |13|14| but planets are initialized with 
their true values (^Cnbody = Ij: 'i'he motion is significantly more 
chaotic, and the peaks are less easily identifiable. Nonetheless, we 
conclude that the primary drivers of Mercury's chaos are Jupiter's 
eccentricity mode and Venus's inclination mode, with Venus's ec- 
centricity mode playing a supporting role. 

narrow range of frequencies in Figure [Tsj Mercury also 
has peaks at |cj| ~ 20'Yyr, due to forcing by Earth and 
Mars. However, these peaks have amphtudes < 10~^, 
and appear to have httle influence on Mercury's chaotic 
motion. Had they been important, one would have ex- 
pected to see their influence in Figure [Tsj whereas all the 
main peaks in that figure have alreadybeen identified. 

5.2. Resonant Angles 

Figure JT6| shows results from a 2 Gyr SWIFT integra- 
tion of the full Solar system (with no reduction of the 
initial e's and i's). The black curve in the top panel is 
Mercury's total \zm\^ which is very nearly equal to its 
total eccentricity (eq. [s]), and illustrates the chaotic be- 
havior of Mercury's orbit. The overplotted green curve is 
the absolute value of Mercury's free zm, which we define 
to be the part of its total zm that comes from the main 
peak in Figure [Tsj i.e., we first take the Fourier transform 
of Zm, then set to zero all frequencies except those sat- 
isfying 4.9'Yyr < uo < 6.5'Yyr, and then take the inverse 
Fourier transform. By plotting the free the short- 
term variations are reduced, and long-term diffusion is 
clearer. The second panel in Figure 16 is the same as the 
top but for Cm ; for the free Cm , we filter out frequencies 
outside of the domain — 6.3'Yyr < a; < — 4.7'Yyr. 

The bottom two panels of Figure [16] show the two four- 
angle combinations involving Mercury that were found 
to undergo libration episodes. The third panel shows the 
angle {wm — ^j) — (^m — ^y), which is the angle that 
has frequency 
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Fig. 16. — Mercury in an N-body simulation of the Solar sys- 
tem: The top panel shows Mercury's \zm\ (approximately its ec- 
centricity) as a black curve, for the duration of a 2 Gyr SWIFT 
simulation. The overplotted green curve is Mercury's free \zm\- 
The second panel shows the same, but for C,m (approximately its 
inclination). The bottom two panels show the two four-angle com- 
binations involving Mercury that were found to undergo libration 
episodes over the course of this simulation. The plotted angles are 
the phases of the free orbital elements (see main text). The angles' 
transitions between libration and circulation are reflected in the 
behavior of of cm and im- 

panel shows the angle associated with the frequency 

cr' {qm - 9v) + {sm - sv) . (56) 




{qm - gj) - {sm - Sv) 



(55) 



i.e., the [1,-1] angle (eq. [43]). (More precisely, we use the 
angles of the free elements; see below.) And the bottom 



Laskar (1992) has shown that the a angle can change 
from libration to circulation. But our finding that the 
can as well is newrj We note that even though both the 
a and a' angles undergo libration episodes, the four- angle 
combination that is the sum of the two, i.e. the angle 
associated with 2gM — gv — 9j does not. That angle was 
found to librate in the A^nbody = 0-75 simulation (Fig. 
13), and we suspect that it will eventually librate in a 
long enough integration of the full Solar system. 

The angles displayed in the bottom panels of Figure 
16 were those of the free elements. For example, for 
wm we first filtered the Fourier transform of zm as de- 
scribed above, and took the phase of the free part of 
zM' This filtering procedure is especially important for 
l]y , because Venus's (y variations are dominated by forc- 
ings due to other modes, including the Mercury-, Earth-, 
and Mars-dominated modes (Laskar|[l990|, whereas we 
wish fly to denote the phase of the Venus-dominated 
mode. Therefore, we first filter the Fourier transform of 
Cy, keeping only frequencies — 7.7'Yyr < u < — 6.3'Yyr, 
and use for the phase of the filtered (y Similarly, 
we filter zm , Cm , i and zj with appropriate windows to 
obtain the other angles of interest. Our method o f filter- 
ing for extracting mode angles differs from that of |Laskar 
( |199Q| ), who extracts mode angles by projecting onto the 
numerically determined nonlinear "proper modes." We 

Since the motion is chaotic, it is possible that a' did not 
librate at all in the simulations of |Laskar| ( |1992] > and [Sussman fc] 
WisdonTI (IT992I . 
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Fig. 17. — Chain of librating angles for the rea l Mercury: This 
is the same N-body simulation as in Figure |16| focusing on the 
time when the [1,-1] angle transitions to its first extended period 
of circulation and then back to libration. Each panel shows the 
angle that has frequency [m,n] = m(gM — 9j) + n(sM — sy), for 
various values of [m,n]. The angles alternately librate, showing 
that the chaos is at least partly caused by the overlap of these 
resonances, as in the /c-model (Fig. [9|. The plotted angles are the 
phases of the free orbital elements. 

have experimented with a number of different methods, 
and also with changing the size of the filter window, and 
the duration of the integration, and found that our filter- 
ing method is simple to implement, is computationally 
efficient, and gives reliable results. 

In addition to the two four- angle combinations of Fig- 
ure [16) one might suspect that there are many more 
higher-order combinations that librate when both of 
those angles simultaneously circulate, as in the Hamil- 
tonian model. In Figure [T7| we zoom into the episode 
when the [1,-1] angle first undergoes an extended period 
of circulation, and plot some higher order combinations 
associated with the frequencies [m, n] = m{gM — 9j) ^ 
Ti{sM — sv)' It can be seen that these angles librate in 
turn. Just as in the Hamiltonian model (compare with 
Fig. [9]). This provides another demonstration that the 
physics of the Hamiltonian model is similar to that of the 
real Mercury. 

6. SUMMARY AND DISCUSSION 

We have shown how secular chaos is driven by the over- 
lap of secular resonances, both for a test particle mod- 
elled with a simplified Hamiltonian, and for the real Mer- 
cury. To linear order, secular frequencies are constant. 
But nonlinearities can shift planets into and out of secu- 
lar resonance with each other, and when two resonances 
overlap, chaos results. 

In §^2][4) we focused on the evolution of a test particle 
in the presence of multiple massive planets. The test par- 
ticle was evolved to leading nonlinear order, and the e's, 
z's, and precession rates of the planets (or more properly 
of the planet modes) were taken to be constant. We first 



considered the simp le case with zero inc linations, as was 
first worked out by Sidlichovsky (1990). In ^ we gen- 
eralized to non-zero inclinations, when the test particle 
comes under the influence of one eccentric and one in- 
clined planet mode. In that case, the particle has two free 
frequencies, its apsidal and nodal frequencies {g and s). 
Each of these is altered by the particle's e and i. There- 
fore each resonance traces out a one-dimensional curve in 
the particle's e-i plane, or equivalently in its Pe-Pi plane. 
A simple way to map out the dynamics is with the "mean 
momentum map" (MMM), whereby the particle's time- 
averaged Pe and Pi are plotted against each other for 
different initial conditions. This shows where the res- 
onances are, how wide they are, and how their overlap 
leads to chaos (Figs. [3][5|. We calculated analytically the 
locations and widths 01 the four strongest resonances — 
the [1,0], [0,1], and [1,±1] — and showed that these agreed 



6]|7| ). Chaos m 



with the numerical MMM results (Figs 
this case emerges from the overlap of resonances of the 
form [m,n] (eq. [43] ), with typically n = m ± 1 and m 
a small integer. Tms may be seen in the MMM, in sur- 
faces of section (Fig. [8|, and also by explicitly tracing the 
chain of librating angles (Fig. [9|. We also examined the 
test particle's trajectories in Fourier space (Figs. |TT|l2). 

In ^ we considered the orbital evolution of Mercury 
in N-body simulations. We showed that despite all the 
simplifications we made in the Hamiltonian models, the 
real Mercury behaved in a qualitatively similar manner. 
In particular: 



• Mercury's chaos is primarily driven by the sy and 
gj modes (i.e. the Venus- and Jupiter-dominated i 
and e modes), although the gy mode also plays a 
role. The nonlinear couplings between those modes 
and Mercury's own free modes (with frequencies 
gM and sm) are primarily responsible for Mercury's 



chaos (Figs. pHTsf . 



• There are a slew of resonant angles that drive Mer- 
cury's chaos. Just as in the Hamiltonian model, 
a chain of resonant angles of the form m{wM — 
wj) + n{Q.M — ^v) show sequential librations, for 



integers [m, n] (Fig. 17), where the angles refer to 
the phases of the free orbital elements. In addi- 
tion, we identified a new four-angle combination, 
{wm — ^v) + (^M — ^v) , that can also undergo 



libration episodes (Fig. 16). 



• Mercury is perched on the threshold of chaos. If 
one reduces the e's and z's of the planets by only 
25%, Mercury's motion becomes nearly regular 
(Fig. 13). This behavior is also apparent in the 
Hamiltonian model (Fig. |5|. We have also per- 
formed a /^nbody = 1-2 simulation, in which the 
planets' initial e's and z's were increased 20% (not 
shown). The result was violent instability, with 
Mercury ejected in ~ 100 Myr. 

Having identified the secular resonances responsible for 
Mercury's chaos, and calculated the widths and locations 
of those resonances, we can calculate the threshold for 
Mercury's chaos. We do that here in an approximate 
way. First, since Mercury's apsidal and nodal frequencies 
differ from gj and sy by ~ 20%, the co-ordinates in the 
Pe-Pi plane where the two resonances overlap are around 
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half that, or pe** ^ Pi** ~ 0.1 (eq. [45]). Second, the 
width of the chaotic overlap zone between the [1,-1] and 
the [1,0] (or [0,1]) is - 2(ejZj) V4(p^^^p.^^)i/8 (eq. 
we include here an extra factor of 2 to account for 



51 
Tlie 

diff eren ce between the half- and full-width, as described 



in ^4.3) Therefore if ej 



^ 3/2 / . 



0.01, then the 



region of chaos will encroach upon the origin of the Pe-Pi 
plane. This explains why Mercury can be chaotic even 
though the eccentricities and inclinations in the Solar 
system are at the level of a few percent. 

The work discussed in this paper can be extended in 
a number of directions. The theory can be extended 
to order-unity eccentricities and inclinations. Although 
that case will be more complicated, we suspect that the 
basic structure will remain, with resonant zones in the e-i 
plane whose overlap leads to chaos. One can also attempt 
to build a theory that includes long-term diffusion and 



massive planets, as well as incorporating MMR's. 

A number of applications also come to mind, such as 
quantifying Mercury's chaotic diffusion and understand- 
ing how it came about that Mercury is perched on the 
threshold of chaos. The latter seems to be a clue for un- 
derstanding how the Solar system arrived at its current 
marginally stable state. It would also be interesting to 
investigate the role of a' (eq. [56]). Is it an unnecessary 
coincidence for Mercury's chaosi 

Our theory can also be applied to Earth and Mars, 
for which librating angles have been identified that are 
similar to those we found for Mercury, i.e. angles of 

the form m(tZ7mars — 'dearth) + ^(^mars — dearth), with 



m,n] = [1,-1] [2,-1] , and [3,-2] (|Laskar|[l992 Suss- 



man & Wisdoml 19921). We also propose that secular 



chaos ca n play a role in shap ing extra-solar planetary 
systems (Wu & Lithwick 2010 ), and hence the theory of 
secular chaos might be applicable to extra-solar planets 
as well. 



APPENDIX 

APPENDIX A: FOURTH ORDER SECULAR HAMILTONIAN 

In this Appendix, we give the expression for the secular Hamiltonian of a test particle perturbed by an external 
planet, where both particle and planet are orbiting a star. The Hamiltonian is expanded to fourth order in the particle's 
eccentricity and inclination, and to leading order in the ratio of semi-major axes. The energy per unit mass of the test 
particle is 

E = -^-^R, (Al) 
2a a' ' ^ ^ 

where Mq is the mass of the star, a and a' are, respectively, the test particle's and planet's semimajor axes, is 
the planet's mass, and R is the disturbing function. We approximate R by only retaining the secular terms up to 
fourth order in e and s = sin(i/2) and second order in a = a /a' (except for the /lo term whose leading contribution 
is 0{a^e^)): 

R^f2e' + fss' + he'e^' + frie's' + eh^' + e^'s') + fss' + Ush^' + hoee^ cos{w - w') 

+ {fiASs' + h^ss'{e^ + e''^) + hQSs'{s^ + s''^)) cos(l] - Q.') + h^e^s^ cos(2ti7 - 29) + /sie^ss' cos(2ti7 - 9' - 9) 

+ hse^s''^ cos(2tu - 29') + f2QS^s''^ cos(21] - 29') , (A2) 

in the notation of the appendix of Murray fc Dermott] (2000). The fi are functions of a that may be expressed as 
sums of Laplace coefficients and their derivatives. We drop terms that are independent of the test particle's orbital 
elements. 

In this paper, we work with a scaled Hamiltonian, H = —2E/^GMQa (eq. Isl), and hence 

(A3) 



H 



^ 7 



3a2 ' 



dropping the Keplerian term in E because it is irrelevant for secular dynamics, and defining 



3 m' 3 /GMo 



1/2 



(A4) 



which is the test particle's secular free precession frequency based on linear theory. The scaled disturbing function 
8i?/(3a^) is a sum of terms that are listed in Tablejl] after expanding the fi to O(a^), and /lo to O(a^). 

APPENDIX B: WIDTH OF THE [1, -1] AND [1, 1] RESONANCES FROM VON ZEIPEL TRANSFORMATION 



We start from Hamiltonian ( 35 ) , which we reproduce here as 

2 2 
P — p- 

H{Pe, qe'^Pi^Qi) = , ' + + ^sPi - ^PePi 



2ejv^cos(7e + 2ij^cosgi , 



(Bl) 



setting 7 = 1, = ^- = — gjt, and Qi = 9_ = 9 — sjt. We solve this Hamiltonian perturbatively, treating cj and 
ij as the small parameters. This is equivalent to expanding in the test particle's forced eccentricity and inclination, 
assumed to be much smaller than the free e and i. We transform to capitalized variables with the von Zeipel generating 
function 

F{Pe,qe;Pi,qi) = Peqe^PiQi + ke{Pe, Pi) sin Qe + ki{Pe, Pi) s'm Qi , (B2) 
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TABLE 1 

^ - Cihi. The Hamiltonian is H = - Cihi. The variables z and C (defined in eqs. 

[8]-[9]) ARE APPROXIMATELY THE COMPLEX ECCENTRICITY AND INCLINATION. TERMS 1-4 ARE Z-ONLY. TERMS 11-17 ARE ( ONLY. TERMS 

21+ ARE MIXED. C.C. DENOTES COMPLEX CONJUGATE. 



Terms in scaled disturbing function ^ 



i 


1 


2 


3 


4 


11 


12 13 


14 


15 




16 


17 


Ci 

hi 


1 


z' + C.C. 


3 

|Z|V|2 


1 

4 


-1 

ICP 


1 1 

C*C^ + c.c. ICPKT 


1 

4 


_ 5 


- C.C. 


_ S 1 
iC^PCc'+c.c. C'C'Vc.c. 


i 


21 


22 


23 




24 


25 


26 




27 


28 




Ci 


-2 




_3 




7 


7 


5 




5 


5 




hi 




\zTh 


\z\'\h 




4 

-C.C. \z'\'^CC + C.C. 


4 

^*2^2 ^ 


- C.C. 


4 


2 

C.C. 2;*^CC^ + C.C. 





where the first two terms generate the identity transformation, and the functions and ki are first order in cj and 
ij] their form will be chosen to "kill" the cosine terms in the Hamiltonian to leading order. The von Zeipel generating 
function transforms variables as follows 



Inserting into the Hamiltonian yields 

p2 _ p2 
Hi^Pe-) Qe-) Pi") Qi) — ^ 

to second order, after setting 



Pe=Pe^ ke COS Qe 
Pi=Pi^ h COS Qi 
Qe=qe -^dp^ke SmQe 

Qi = qi -^dp.kisinqi . 



APe + AgPi - 2PePi - keki COs{Qe - Qi) - keh COs{Qe + Qi) 



ke — 

ki = 



2ejVPe 



A - Pe/2 - 2Pi 
-2ijVPi 

As + Pi/2 - 2Pe 



(B3) 
(B4) 
(B5) 
(B6) 



(B7) 

(B8) 
(B9) 



to eliminate the first order terms. The two cosine terms in this Hamiltonian are the [1, -1] and [1, 1] resonances. 



respectively (see eqs. 41 - 42 and following). Note that we have dropped second order terms in the Hamiltonian 



that are proportional to cos^Qe, sin^ Qe, cos^ Qi, and sin^ Qi, because these have little effect on the [1, 1] and [1, -1] 
resonances; but they generate new frequency components which are important for higher order resonances. 

To leading order, ke is twice the product of the free eccentricity (v^) with the forced eccentricity, where the forced 
eccentricity differs from the linear expression (ej/A; see e^ [18]) by the terms Pe/2 + 2Pi in the denominator, which 
arise from the nonlinear shift of the frequency (eqs. 42 )^imilarly ki is twice the product of the free and forced 
inclinations. Hence the strengths of the [1, -1] and [rji] resonances are proportional to the products of the free and 
forced eccentricities and inclinations, as argued qualitatively in Q 

To determine the width of the [1, -1] resonance, we drop the last cosine term in the above Hamiltonian. Since 
P+ = Pg + Pi is an integral of the motion, we may re-write the Hamiltonian as i^(Pe, Q-) = 2(Pe — P*)^ — keki cos Q_, 
dropping a constant and defining Q_ = — Qi and P* = (5/8)P+ — (1/4)(A^ — A^). Therefore the half- width of the 
resonance is 

SPe = ^/\Kk~\ . (BIO) 

We take the amplitude of the cosine term to be fixed at its value at resonance center (e.g., Chirikov 1979| ). Since 
Pg + Pi is constant, the half-width in Pi is the same, 6 Pi = ^/\k^ki\. 
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